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ABSTRACT 


We present a detailed investigation of the flaring activity observed from a BL 
Lac object, S5 0716+714 , during its brightest ever optical state in the second 
half of January 2015. Observed almost simultaneously in the optical, X-rays 
and y-rays, a significant change in the degree of optical polarization (PD) and a 
swing in the position angle (PA) of polarization were recorded. A detection in 
the TeV (VHE) was also reported by the MAGIC consortium during this flaring 
episode. Two prominent sub-flares, peaking about 5-days apart, were seen in 
almost all the energy bands. The multi-wavelength light-curves, spectral energy 
distribution (SED) and polarization are modeled using the time-dependent code 
developed by Zhang et al. (2014). This model assumes a straight jet threaded 
by large scale helical magnetic fields taking into account the light travel time 
effects, incorporating synchrotron flux and polarization in 3D geometry. The 
rapid variation in PD and rotation in PA are most likely due to re-connections 
happening in the emission region in the jet, as suggested by the change in the 
ratio of toroidal to poloidal components of magnetic field during quiescent and 
flaring states. 
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Introduction 


Blazars are an extreme subclass of active galactic nuclei (AGN) known to posses an 
extremely collimated relativistic jet, perpendicular to the plane _of the accre tion disk and 
seen at very small angles (<15°) to our line of sight (LOS) (Urry & Padovani 1995, 120001 ). 


Such close alignment of the jet to the LOS leads to the relativistic boosting of the jet emis¬ 
sion, which dominates the blazar emission. The spectral energy distribution (SED) of a 
blazar typically exhibits bi-modality with two broad bumps: one peaks between sub-nun 
to UV/X-rays, and the other peaks somewhere at the MeV-GeV energies. The low en- 
ergy part o f SED is well established to be due to synchrotron process in the relativistic jet 
flUrrv fe Mushotzkv 198 2). with other non-jet components like disk, torus, BLR etc. also be¬ 


ing s i gnificant contributors in several sources albeit in narrow energy range (jGhisellini et al. 


20091: iNalewaiko et al.ll20l4 iKushwaha et al.l I2014J) . The high energy bump, on the other 
hand, is poorly understood and according to one approach, this component is thought to 
arise from inverse Compton scattering of low energy seed photons by the highly energetic 
leptons (e + /e _ ) in the jet. If the synchrotron photons, originating in the same population of 
high energy lepto ns are upscattered to high energies, the process is termed as synchro tron self 
Compton (SSC) ( Bloom fe Marscher 1996 : Sokolov et al.l 2004 : Ghisellini et al. 1985). In the 
cases where the seed photons are external to the jet, namely, disk, torus, BLR and/o r some¬ 


times even CMBR, then the process is known a s external Comptonization (IDermer et al. 
1992 : Sikora et al. 1994 . 2009 : Agudo et al. 2011). 


The emission from blazars shows enormous amount of variability at almost all the fre¬ 
quencies. Since the central engine is not resolvable by any existing facility, the variability 
provides a useful tool to diagnose the physical mechanisms responsible for the emission, 
thanks to the availability of quasi-simultaneous data at various energies from space-based 
(e.g. Fermi, Swift etc. ) and many ground based observatories. Optical polarization obser¬ 
vations can sufficiently constrain many of the jet properties e.g., the strength and nature of 
the magnetic field, the geometry of the jet, and the physical processes etc., provided these 
observations are properly supplemented by the variability information at other energies. The 
high energy y-ray flares are mostly followed by the activities at lower energies with a few 
exceptions of orphan flares. The observability of the emission at a particular energy from 
a blazar jet may be dependent on the opacity of the emission region at that wavelength. 
Therefore, in this scenario, the location of the dissipation region becomes an important 
factor behind the nature of variability seen across the electromagnetic spectrum. In other 
words, a multi-wavelength variability study of blazars can also prov i de us information a bout 
the location of the emission region in the jet (IMarscher et al.ll2Q10l: lOrienti et al.ll2013l ). 


In the literature, some cases of rotation or swing in PA have been reported during high 
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energy flares in blazars flAbdo et al.ll2010t lMarscherll20l4 IZhang et al.lI2Q14J) . Being depen¬ 
dent only on the orientation of shock and magnetic field threading it, PA provides an unique 
tool to understand the acceleration mechanisms and the behavior of the shocked plasma. 
Several models have been proposed to understand the PA swings, however, currently only 
two of them are capable of reproducing the multi-wavelength light-curves, time-dependent 
SEDs, and multi-frequency polarization, simultaneously. The first mod el, known as the Tur¬ 
bulent, Extreme Multi-zone Model (TEMZ), suggested by IMarschen (20141). assumes that 
the emission region comprises of turbulent chaotic magnetic field with a large number of 
small regions of ordered small scale magnetic fields called ‘cells’. A proper accounting of 
time dependent contributions of individual cells to the total emission is done to simulate the 
observables. The seco nd a pproa ch termed as the Helical Magnetic Field Model (HMFM), 
suggested by Zhang et ah (12014]) assumes a large scale ordered helical magnetic field with a 
more rigorous and proper accounting of the Light Travel Time Effects (LTTEs) for individual 
zones to recreate the SEDs, light-curves and polarization. 


The blazar S5 0716+714 , at a red-shift, z of 0.31 ( Nilsson et al.l 20081 ) is one of the 
brightest BL Lac ob jects that is highly variable from radio to y-rays with a very high 
duty cycle (IWagner et al.l 1 9961 ). In the SED based classification sch eme, S5_0716+714 is 
sub-classified as an intermediate energy peaked BL Lac object (IBL) (Giommi et al. 1999), 
and is also confirmed based on the con cave shape of its X-ray (0.1-10.0 keV) spectrum 
(iFoschini et al.l 120061 : iFerrero et al.l |2006|) . The concave shape of the spectrum in X-rays is 
indicative of the presence of a tail from the synchrotron emission (falling tail) and a flatter 
part from the inverse Compton (IC) spectrum (rising tail). This object has been detected 
in the MeV energy range several times at different flux levels by the Energetic Gamma Ray 
Experime nt Telescope (EGRET) detector on board the Compton Gamma-ray Observatory 
( CGRO ) ( Hartman et al. 1999 : Nandikotkur et al. 2007). In 2008, AstroRivelatore Gamma 
a Immagini Leggero {AGILE) reported the detection of variable y-ray flux with a peak flux 
density above the maximum reportec 
in the Fermi-LAT bright source list 


by EGRET (IChen et al.l 120081) . S5 0716+714 is also 
Abdo et al. 2009lh MAGIC first detected this source 


in VHE y-rays during 2007 November (A>o. 4 TeV 
much larger flux (F >0ATe y_ 


0.8 x 10 11 erg cm 2 s x ) and later with 
7.5 x_10 -11 erg cm ~ 2 s -1 ) at 5.8 a, in 2008, April. Similar trend 


was also seen in optical (lAnderhub et al.ll2009l) while source was in historically high state in 
X-rays ( Ciommi et ab 20081) in 2008, April. A concurrent rapid rotation of the PA w as also 
observed just after the maximum in optical flux had reached (lLarionov et al.l 20081) . This 
seems to sup port the indication seen in the previous MAGIC observations for other BL Lac 
objects (lAlbert et al.l 12006 . 2007). that there is a connection between the optical high states 
and the VHE y -ray high states. 


The recent optical monitoring of S5 0716+714 shows a consistent rising trend in the R 
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band flux over the past few months. The older observations reveal that the source had gone 
to a very faint state (14.9 mag in R-Band) around 2013, December (MJD 56650), very close 
to the fain test s t ate o f this source reported till date (private communications). Very recently, 


Carrasco et ahj ([2015) had reported a high state in IR on MJD 57033.3 (2015, January 11) 


which was around 2.5 magnitude brighter than the previously observed flux on MJD 57021 
(2014, December 29). Following these observations, there were several teleg rams r eport ing 
a further brightening of S5 0716+71 f in different optical and IR bands 


2015 


Bachev et al.ll2015t iBachev fc Strigachevll2015l: ISpiridonova et al. 12015 


Arkharov et ah 


Chandra et ah 


2015). Soon after these reports of optical/IR flares, Mirzoyan (120151 ) reported a variable 


VHE detection above 150 GeV. This VHE detection, quasi-simultaneous with a high flux 
state in the optical, seems to be similar to that seen in the 2008 event. In the present 
work, we have investigated this event in the framework of time-dependent modeling of the 
observables, namely, light-curves, SEDs and multi-frequency polarization, to understand the 
role of the magnetic field in the blazar jet. Our study is mainly focused on the simulation 
of the part of the outburst where simultaneous data at all energies (optical to y-rays) are 
available. This paper is organized as follows: §2 provides the details of data resources and 
the analysis methodology adopted. We present the multi-wavelength light-curves and SEDs 
along with the modeling of these in §3, followed by the discussion of the results obtained, in 
§4. Finally, in §5 we present a summary of our work. 


2. Observations &; Data Analysis 


We used data from Large Area Telescope (LAT), on-board Fermi spacecraft (Atwood et al. 


2009), for y-ray counterp arts of the optical fl aring event. The Swift/X .RT jBurrows et al. 


20051 ) and Swift /UVOT (iRoming et al. 2005 ) data are analyzed for X-rays and UV light- 
curves and SEDs. The optical monitoring data from 1.2 m & 0.5 m telescopes at Mt. Abu 
Infra-Red Observatory (MIRO) are also analyzed. A few other observations, reported in 
“The Astronomers Telegrams”, are also used in the present study. The corresponding values 
are corrected for galactic extinction and reddening. The publicly available spectropolari- 
metric observations from Steward Observatory, Arizona, are also used. In the following, we 
summarize the details of the observations and data analysis techniques used for various data 
sets. 
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2.1. Fermi-LAT 


LAT normally works in all sky scanning mode, covering the whole sky every 3 hours. The 
scanning mode along with its large held of view (FOV) provide almost 30 min of monitoring 
of each source during the course of each scan. The broad energy coverage of LAT (A E ~ 
0.02-300 GeV) makes this facility ideal for studying high energy astronomical events. 


Two months of Fermi-LAT P7REP data from 2015, January 01 to 2015, March 01 were 
analyzed using the Fermi ScienceTools v9r33p0. In the analysis, only SOURCE class events 
with energies between 100 MeV - 200 GeV, defined under the Instrument Response Function 
(IRF) P7REP_SOURCE_V15 , have been used from a 15° region of interest (ROI) centered 
at the location of S5 0716+714 (R.A.=110.473°, DEC.=+71.343°). We have not included 
the low energy and high energy ends of LAT data to avoid the possible contamination with 
artificial counts cau sed by the poor response of the detector in this range. All the sources 
from 2FGL catalog ( Nolan et ah 2012) within the ROI, plus an annular radius of 10° around 
it, were modeled along with the standard diffuse templates of Galactic (glLiemjv05jrevl.fit) 
and isotropic background (is o sources 05-rev 1.txt) . For generating the XML model hie, we 
have made use of the contributory python package, make2FGLxml.pj0. The criterion for 
using a proper model for generating light-curves a nd SEDs is similar to the one adopted by 


Chandra et al.l (120141 1 and Kushwaha et al.l (120141 1 . At first, all the sources within ROI+10 C 


were considered for unbinned likelihood analysis. The point sources with TS values less than 
0, obtained from likelihood analysis, were removed from our input model. The likelihood 
analysis procedure was repeated until it converged without any source with TS < 0. The 
daily fluxes for the lightcurve were then extracted using the best model parameters where 
both spectral indices and normalization were kept free. However, for generating the SEDs, 
spectral indices in the model were kept fixed at -1.6 ± 0.01 as obtained from the likelihood 
analysis of the complete data set. 


2.2. Swift Data Analysis 

The 0.3-10.0 keV X-ray and UV/optical archival data from XRT and UVOT instru¬ 
ments, on-board the Swift X-ray space-borne observatory ( Gehrels et al . 2005), were ana¬ 
lyzed for the present work. The standard procedures prescribed by the instrument teams 
were followed step by step to generate the science products. The recently updated version 
of the calibration database (CALDB) along with heasoft v6.16 were used for our analysis. A 


1 http://fermi.gsfc.nasa.gov/ssc/data/analysis/user/ 
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number of pointings made by Swift during 2015, January 19-31 were analyzed. The following 
is a quick summary of the analysis procedures used for XRT and UVOT data. 


2.3. XRT Data 


The level 2 cleaned event hies were generated using the standard xrtpipeline tool with 
the default parameter settings and by following the prescriptions of the instrument team. 
The source and background lightcurve and spectra were generated with appropriate region 
and grade filtering using the xselect tool. In this case, the source spectrum was extracted for 
a circular region of 47" radius around the source location while four source-free regions in the 
neighborhood of the target, each with a 100" radius, served for the background spectrum. 
The count rate in these observations exceeded the recommended pile-up free count rate for 
PC mode (0.5 counts/sec). Therefore, all the event hies were investigated for pile-up and a 
proper procedure for pile-up correction as suggested by other researcher^] were used wherever 
needed. The ancillary response matrix was generated using the task xrtmkarf followed by 
the xrtcentroid task. The response matrix hie provided with the CALDB distribution was 
used for further analysis. 


Spectral fitting was done in the energy band between 0.3 to 10.0 keV using the XSPEC 
(version 12.8.2) package distributed with heasoft package. A simple power law along with 
the Galactic absorption gives the best ht for almost all the observations of interest. The 
model parame ter, Nh i .e., the interstellar column density, was kept hxed at a value of 3.06 
x 10 20 crrr 2 3 ( Kalberla et al. 2005). This value jvas estimated using the web-based tool 


developed by a group at the University of Bonijf). The normalization and spectral index 
of the power law were the free parameters for the spectral htting. Table 1 summarizes the 
values of various parameters obtained from the spectral htting for different observation IDs. 
The 0.3-10.0 keV huxes thus obtained were used for constructing the light-curves. 

The Galactic absorption co rrected X-ray energy spectrum was constructed using the 
procedure adopted in Chandra et al. (2014) with that instead of using the default binning, 
all spectra hies were binned according to a hxed input hie (describing the details of channel 
binning) using the grppha tool. The multiple X-ray spectra, if any, for a particular SED 
segment, were merged keeping in mind the possible spectral variations. 


2 htt p: / /www. swift. ac. uk / analysis/xrt/pileup. php 

3 https://www.astro.mii-bonn.de/hisurvey/profile/index.php 
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2.4. UVOT Data 


UVOT snapshots with all the six available filters, V (5468 A), B (4392 A), U (3465 A), 
UVW1 (2600 A), UVM2 (2246 A), and UVW2 (1928 A) for all the ObsIDs were integrated 
with the uvotimsum task and analyzed using the uvotsource task, with a source region of 5", 
while the background was extracted from an annular region centered on the source location 
with external and internal radii of 40" and 7", respectively. The fluxes thus obtained, were 
corrected for galactic extinction using a too l developed for R platform. This interactive 
tool adopts the model described in Carde lli et al. (119891) . This tool needs E(B-V) value as 
input, which was esti m ated us ing the web-based calculator by NASA/IPAC Infrared Science 
Archive^ ( Schlaflv V Finkbeiner 2011 ). The value of this parameter in the direction of S5 
0716+714 is 0.026 ± 0.001. The extinction corrected fluxes obtained from the aforementioned 
procedures were then used for extracting the light-curves and SEDs. 


2.5. Optical Photometry & Polarization Data 


The photometric monitoring of S5 0716+714 was performed using iKon ANDOR CCD 
Camera (2048 x 2048) as a backend instrument at the f/13 Cassegrain focus of the 1.2m 
optical telescope of MIRO. The CCD in this camera is cooled to —80° with thermoelectric 
(TE) cooling to keep the dark current very low. The camera is attached to a coupling unit 
consisting of a filter wheel with 12 slots, equipped with 10 optical filters (Standard John¬ 
son/Cousins UBVRI broadband filters + u,g,r Sloan filters + two narrow band filters) just 
on top of the camera. One of the two additional slots is kept blocked for grabbing bias/dark 
frames, whereas the other is kept open for white light monitoring of extremely faint sources. 
The bias and sky flats were taken on daily basis for performing the pre-photometric image 
processing. The source was also monitored using the Automated Telescope for Variability 
Studies (ATVS) at MIRO mounted with TE cooled (~ —80°) ANDOR iXon EMCCD camera 
(1024 x 1024) as the backend instrument. This camera is coupled with a filter wheel with 
the Standard Johnson/Cousins UBVRI filters. The twilight sky flats and bias frames were 
also captured on daily basis. The observing strategy for both the facilities was to capture 
a few images in all the bands (UBVRI) and then monitor for a longer time in one filter, 
say R band. The apparent magnitudes of the source and the known field stars (photometric 
standards) for individual exposures were derived using th e sta ndar d aper ture photometry 
technique preceded by pre-processing of images (iChandra et ah 2 0111 2014). The magnitude 
was then converted into flux using simple conversion factors and zero point flux (ZPF) for 


4 http: //irsa.ipac.caltech.edu/applications/DUST / 














different bands (IB esse ll 1979). Apart from photometric monitoring at MIRO, the observed 
fluxes from various Astronomers Telegrams were also used for comparison and completeness. 
The R band fluxes obtained from the Steward observatory database are also used followed 
by extinction correction. 

Optical polarimetric data used in this study were taken from the observations performed 
as a part of the Fermi Support observing Program at Steward Observatory, Arizona, USA 
flSmith et al.l 20091 ). The PD and PA values provided at Steward observatory data base are 
already calibrated and hence can be readily used. We have corrected the PA values for 180° 
angle ambiguity, in order to make the rotation clearly visible (Fig. |T](f )). 


3. Results and Interpretation 

The multi-wavelength light-curves derived using observational data from various re¬ 
sources, mentioned in §2 are shown in Fig. [0 The SED for the duration of MJD 57045.5- 

57047.5 was also generated using data from the aforementioned facilities. In the following, 
we discuss the results and their interpretations. 


3.1. Multi-wavelength light-curves 

The Figs. 1(a) & 1(b) show the y-ray and X-ray fluxes observed by Fermi & Swift , 
respectively. The Fig. 1(c) shows the UV/optical light-curves derived from UVOT obser¬ 
vations while Fig. 1(d) shows R band magnitudes from various resources. The Figs. 1(e) 
& 1(f) are the PD and PA observations from the Steward observatory. The figures show 
significant variability in all the energy bands (y-rays to UV/optical) as well as in the optical 
polarization. As evident from Fig. 1(a), y-ray show a clear trend of multiple ups and downs 
during the span of MJD 57034 to 57055. A consistent rise in flux is seen between MJD 

57034.5 to 57039.5, reaching to a flux level of 0.75 x 10~ 6 ph cm~ 2 s _1 , which is about 4 
times the average y-ray flux of this source. A rapid variation in flux is seen during MJD 

57040.5 to 57042.6 reaching to a flux level of 0.83 x 10 -6 ph cm -2 s -1 on MJD 57041.5, i.e., 
within one day, the flux increased and then decreased to 0.35 x 10~ 6 ph cm -2 s^ 1 the next 
day. The flux again started rising slowly after MJD 57042.5 and reached upto 1.1 x 10~ 6 
ph cm -2 s _1 , the highest y-ray flux level reported for this source. However, a fast decrease 
in flux soon after this flare is also seen. We, therefore, notice two major y-ray sub-flares 
associated with the January 2015 major flare in S5 0716+714 ■ The Swift X-ray light-curve 
in 0.3-10.0 keV band (Fig. 1(b)) exhibits a consistently rising trend with a peak flux of 
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2.7 ±0.3 x 10 11 erg cm 2 s 1 , on MJD 57047.2, which is comparable to th e pea k X-ray flux 


seen during 2007 October-November flare of the S5 0716+714 f Oiommi et ah 2008). The 
poor coverage of the XRT observations restricts us from making conclusions about multiple 
sub-flares in X-rays. However, the ups and downs in the individual data points provide a 
glimpse of small scale variations. 

The R band optical light-curve (Fig 1(d)), shows two distinct, well separated sub-flares. 
Initially, the R band flux slowly rises and reaches to 11.6 mag from 12.41 mag during MJD 
57035 - 57041.3 (rate ~ 0.13 mag /day). The peak flux here corresponds to brightest ever 
state of the source reported so far fjChandra et al. 20111). Just after this peak, the flux 
decreases to 12.22 magnitude between MJD 57040.5 to 57044.4 (0.14 mag/days). Here, we 
have ignored the variations in subsequent nights. Specifically, the source had undergone 
very fast decrease in flux (AM ~ 0.5 mag) during MJD 57040 to 57041.8 and then again 
went to 11.9 mag on MJD 57042.4 before falling to 12.2 mag on MJD 57044.2. Soon after, 
the R band flux again starts rising at a slightly faster rate (0.33 mag/day). The peak flux 
corresponding to the second bump is almost equal to that of the first one. However, the 
flux value of the second bump remains constant for quite some time (~ 2 days). After this, 
the flux gradually decreases to 12.68 mag within MJD 57046 to MJD 57054.5 (rate ~ 0.1 
mag/day). The later monitoring suggests even fainter state of this source. Note that the 
plotted flux values represent fluxes averaged over few hours of monitoring. The fluxes in the 
other optical bands, namely, V, B and I, also show the similar behavior, but are omitted 
in the light-curves for the sake of clarity of Fig. [0 Fig. [21 depicts the correlation between 
fluxes in different optical bands and the flux-dependence of the B-V color. The fluxes in 
various optical bands are differently correlated, which is clearly seen from the slopes of the 
best linear fits to these curves (Fig. [2]). The second block of Fig. [2] represents the standard 
‘bluer when brighter (BWB)’ trend seen during a typical flare in blazars. The Swift UVOT 
observations (Fig. 1(c)) also indicate the existence of two humps in the light-curve. However, 
the poor coverage by Swift pointing provides an incomplete picture of the flux variations. 


The variations in the optical PD and PA are presented in the Figs. 1(e) & 1(f), respec¬ 
tively. The PD shows various episodes of rapid variation with the trends completely unrelated 
to the total flux variations. More specifically, at the very beginning of our available dataset 
i.e., around MJD 57042, when S5 0716+714 had already passed through the first optical 
sub-flare, the PD was very high (10.67 ± 0.02 %), which then decreased to 4.0 ± 0.02 % by 
MJD 57045.3. The next two observations reveal a 5.4 % change in PD between two epochs in 
the same night differing by 2.4 hours (MJD 57045.3-57045.4). The next episode of PD is even 
more dramatic and nicely covered as the PD decreases by 4.68 % within ~ 7.2 hours (8.88 
- 4.20 % between MJD 57046.2 to MJD 57046.5). During the next three segments, namely, 
MJD 57047.2-57047.5, MJD 57048.1-57048.2 and MJD 57050.3-57051.5, independent rising 
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trends of 4.4 % (6.38-10.79 %), 2.96 % (8.0-10.96 %) and 8.01 % (0.45-8.46 %), respectively, 
are seen. The typical error in these observations of PD is 0.02 %. A straight line was fitted to 
the individual segments of PD variations using least square fitting algorithm, which clearly 
indicates very different slopes for all of them (LI to L5 in left panel, Fig. El)- It is very 
difficult to associate rapid fluctuations in PD to the double-humped shaped variation in the 
total flux light-curve. It could perhaps be due to two emission components contributing to 
the total flux; one nearly unpolarized and other polarized; both varying with time. The 
observed PA profile also supports the same argument, as discussed below. Initially, during 
the first segment of the PD variations (MJD 57042.3-57045.3), the PA is mildly variable 
around 100° with a slight change from 78.6° to 100.2°. This implies the decay of a previous 
polarized component and the emergence of a new one; both with similar PA. Afterwards, 
the new polarized component starts to evolve. During MJD 57045 to 57047.5, the PA con¬ 
sistently increases by ~ 164°, from 93.0° to 256.9°. The later observations (MJD 57048.1 
to 57051.5) suggest a consistent decrease in PA from 256.9° to 96°. The rising and falling 
parts of the PA variations were fitted with an exponential function [f(x) r\j 4 e -a(*-x 0 )], tQ 
get an estimate of the temporal profile (right panel, Fig. [3]). The index of the exponential 
for the rising part is —0.45 ± 0.01 day _1 whereas that for the falling part is 0.29 ± 0.006 
day _1 . This shows that the second PA swing is slower than the first one (right panel, Fig. 
[3]). These make the case for two rotations, of about 180° each, during the course of the 2nd 
optical sub-flare, indicating the sub-flare to be strongly related to a significant change in the 
magnetic field. 


Mirzoyan (2015) reported MAGIC detection of the variable VHE y-ray flux in the range 
of 4 x lO^ 11 to 7 x 10 -11 ph cm -2 s -1 above 150 GeV, during MJD 57044-57048, which is 
simultaneous with the X-ray and optical high states of the S5 0716+714 during the second 
sub-flare (Fig. [[]). It appears that this variable VHE detection is correlated with the activity 
in the optical and the X-ra ys, s imilar to the 2008 detection of VHE y-ray emission of S5 
0716+714 (lAnderhub et al.l 20091 ). In the following, we have described the modeling of the 
observables derived for the duration of MJD 57045.5-57047.5. 


3.2. Modeling of Light-curves, Optical Polarization SED 


In general, there are two mechanisms, namely, shocks and magnetic reconnection, that 
ma y re sult in &_ flaring activity in a blazar jet. We have used the HMFM model proposed 
by f Zhang et al. 2014 1 to fit the first PA rotation, from MJD 57045.5 to 57047.5 (the rising 
part). This model assumes an axis-symmetric cylindrical geometry for the emission region, 
which is further evenly subdivided into zones in radial and longitudinal directions (see Fig. 
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Zhang et alj 2014). The relativistic plasma in the emission region, pervaded by a helical 


magnetic held with a possible additional turbulent component, moves on a straight trajectory 
along the jet and encounters a hat stationary disturbance. The disturbance will temporarily 
modify the physical conditions (magnetic held, particle distribution etc.) at its location 
inside the emission region. The initial state in the emission region is regained once it moves 
out of the shocked region. The non-thermal particle evolution, radiation and polarization 
signatures of our model are realized by the M onte-Carlo/Fokker-Planck (MCFP) radiation 
transfer code developed by Chen et ah f 2012h and the 3D polarization dependent radiation 


transfer code developed by IZhang et al . (2014). In this model, the particle distribution 


in individual zones are evolved using a locally isotropic Fokker-Planck equation, and the 
polarization dependent emission from each zone is properly traced to account for all the 
LTTEs. Even if the disturbance in our model is hat, the haring region, which consists of the 
zones affected by the disturbance whose emissions arrive to the observer at the same time, 
will be distorted into an elliptical sh ape. Illustrations of the model setup and the LTTEs 
can be found in Fig. [4] and also in Zh ang et al. (120151 ). 


Table 2 lists some key parameters from our modeling, estimated in the co-moving frame 
of the emission region. Due to the relativistic aberration, even though we are observing 
these objects at very close alignment to the relativistic jet in the observer’s frame (typically, 
0* bs ~ 1/r, where F is the bulk Lorentz factor of the outhow along the jet), the angle 6hbs 
between LOS and the jet axis in the comoving frame is likely to be much larger. Specifically, 
if 6* bs = 1/r, then 0 o b s = 90°. In our fitting, we choose the LOS in the comoving frame 
at $obs = 90°. As is mentioned above, the PD variations are hardly linked to the flux 
variations; therefore, an unpolarized turbulent contribution to the total flux is likely to be 
present. Therefore, we have fitted the averaged SED for the period, which contains a helical 
component and a constant turbulent contribution, then fit the polarization signatures based 
on the derived parameters. Also the PA profile is fitted prior to the PD profile, as the former 
is less affected by the turbulence. The best fit model and data are displayed in Fig. EjJ 


We assumed a leptonic origin for the SED, so that the low-energy component is domi¬ 
nated by synchrotron, while the high-energy component consists of both SSC and EC con¬ 
tributions. Due to the Monte-Carlo photon tracing, small numerical errors are present, but 
the overall fitting is very close to the data. The time-dependent polarization profiles are 
interpreted as alterations in the magnetic topology initiated by the disturbance. The ori¬ 
gin of the disturbance can be either shock or shock-initiated magnetic reconnection. Before 
the emission region encounters the disturbance, the entire region contributes uniformly. Al¬ 
though the helical magnetic held has comparable poloidal and toroidal contributions, due to 
the axis-symmetry of the emission region and the LOS orientation, the projected poloidal 
contribution onto the plane of sky is stronger, resulting in a poloidal-dominating polarized 
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flux. When the disturbance moves in, it will alter the local magnetic field to be toroidal- 
dominating. In addition, it will strengthen the emission in this region by amplifying the 
local magnetic field strength (shock) or injecting additional non-thermal electrons (recon¬ 
nection). Due to the LTTEs, only the near side of the flaring region is observed at first, 
so that the PD gradually drops as the initial poloidal dominance is balanced by the flaring 
toroidal contribution, and the PA rotates towards toroidal component. At the middle of 
the event, the flaring region will extend across the emission region, therefore it will receive 
equal contributions from both the near and the far sides (Fig. 0]), as seen from the ob¬ 
server, mimicking the initial axis-symmetry, but dominated by toroidal contribution. Since 
the length of the cylindrical emission region is longer than its diameter, this flaring region 
will stay for some time, creating a “step phase” in both PD and PA profiles, as proposed in 


Zhang et al. (20151). After that, the disturbance will completely leave the emission region, 


and the flaring region moves to the far side, so that the PD reverts back to its initial state 
in a time-symmetric pattern. The PA instead completes a 180° rotation to the initial state 
(notice the 180° ambiguity), as the projected toroidal component on the near side is opposite 
to that on the far side. 


We have further investigated a few issues in the polarization fitting. We notice that 
although the decreasing part of the PD profile is well reproduced, the rising part around 
MJD 54047 is a bit off. We remind the readers that we have applied a constant unpolarized 
turbulence in the fitting, which in reality, may actually decrease with time, consequently 
PD may increase. The PA, however, is generally unaffected. Moreover, the “step phase” is 
indeed necessary and can put a strong constraint on the ratio of the length and diameter of 
the emission region. We can see in the data that both PD and PA tend to converge at a 
stable value during the middle of the event. Most importantly, the slopes of the PA rotation 
before and after the “step phase” set a stringent constraint on the ratio of the flaring toroidal 
dominance to the initial poloidal dominance. This ratio not only leads to a definite excess 
of toroidal contribution during the “step phase”, which is shown in the PD profile, but also 
require the helical component to undergo an increment in the flux. This is consistent with 
the flare activity observed at X-rays and y-rays. During the PA rotation, the total optical 
flux is almost stabilized to its highest observed flux and hardly shows any pattern. Therefore, 
to compromise the flare in the helical component, the turbulent contribution has to decrease, 
further confirming our previous assumption. 
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4. Discussion 


The observations of simultaneous variations in multi-wavelength bands shown in Fig. 
l(a)-l(f) suggest that the similar activity is driving the emission over a broad range of 
spectrum during the brightness phase of S5 0716+714 ■ The fact that the flare is seen almost 
simultaneously over a broad spectrum, further suggests that the emission region is likely to 
be co-spatial. Additionally, it indicates that the region of maximum dissipation should be 
transparent to the emissions at all energies. The variable VHE detection that followed the 
activity in optical along with the historically bright state in X-rays and y-rays, coinciding 
with variable PD and rotation of PA, imply that the same electron population is responsible 
for this major activity. In short, in our model, these observations support a leptonic origin 
for the high energy emission during the 2015 activity of S5 0716+714 ■ A careful look at the 
Fig. CDa) " IHf) reveals occurrence of two major sub-flares, in all the bands, super-imposed 
with small amplitude fluctuation. Any analysis of the time lag between different wavebands 
is, however, difficult to perform due to limited coverage and their resolution in the present 
data. The BWB trend and the different slopes in flux-flux correlation graph (Fig. [2D indicate 
that this flare may have been caused by shock acceleration activities in the jet rather than 
being triggered by the involvement of geometry dependent effects. The very rapid (hourly) 
and prominent variations in PD reflect the crucial role played by the magnetic field during 
this event. The individual segments of PD variations may be due to the fluctuations in the 
shocked region resulting in changes in the magnetic field in the compressed region. 


The variations in PA observed during the second optical sub-flare is the most important 
feature of this state. There is a clear indication of rotation in PA, almost coincident (within a 
day) with the optical, X-ray and y-ray flares. This particul ar scenario is well observed in some 


of the blazars namely 3C279, Mrk421, PKS 1510-089 etc. flAbdo et al.ll2010t iMarscher et al. 


2010l ; lMarscherll2014j) . Several mechanisms have been proposed to interpret the PA rota tions, 


such as an emissio n region moving alon g a curved trajectory (jVillata fe Raiteril Il999h or a 
bendi ng in the jet (iMarscher et al.lll99ll) . or streamlines following the helical magnetic field 
lines I Marscher et al.l 2008 1. or stochastic activation of individual zones in a turbulent shock 
(TEMZ, model). The bending jet model involves an asymmetric jet structure and thus 
requires a lot of freedom in the parameter space. Streamlines along a helical magnetic 
field imply that the PA rotation should be preferably in the same direction, while the two 
sequential PA rotations are observed in the opposite directions. Additionally, both models 
cannot naturally explain the simultaneous flaring activities. The TEMZ model, on the other 
hand, can hardly produce the systematic, apparently time-symmetric polarization profiles 
due to its stochastic nature. Nevertheless, our HMFM model applies a simple axisymmetric 
geometry for the emission region and the disturbance, taking into account all the LTTEs, 
naturally explaining the simultaneous flaring activities and the apparently time-symmetric 
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polarization profiles. In addition, based on the behavior of the disturbance, an immediately 
following PA rotation event in the opposite direction is possible. However, since the second 
PA rotation does not have sufficient data coverage, we are unable to constrain the model 
parameters. This part is omitted in Fig. 0 


Our best fit parameters do not allow us to distinguish between the shock scenario and the 
shock-initiated reconnection scenario. However, we still prefer the reconnection mechanism 
because our model predicts higher flare amplitude at high energies than in the optical, as is 
seen here. In the reconnection, the increase of the non-thermal electron density is due to the 
dissipation of magnetic energy, which weakens the synchrotron flare. The shock scenario, 
however, enhances both the non-thermal density and the magnetic held strength, leading to 
a stronger synchrotron hare. The aforementioned hare and PA r ota tion event i s very s imilar 
to the other hare + PA rotation event s seen in blazars flAbdo et al.l l2010l: iMarscher et al. 


2010l ; iMarscherl 12014]; IZhang et al.ll2015l) . After these events, hux in all bands suffers a large 


decrease and becomes less active. This indicates that a severe energy dissipation occurs 
during the hare+PA rotation. The PA rotation implies a strong alteration in the magnetic 
held, which again provides another piece of evidence for the reconnection. 


5. Summary 

S5 0716+714 exhibited a multi-wavelength outburst in January 2015, with two well- 
resolved sub-hares in gamma-rays and optical. The event was accompanied by rapid PD 
variations and a systematic PA swing, and followed by a variable VHE detection (F >150 ceV 
= 4-7 x lO^ 11 ph cm -2 s -1 ) by the MAGIC group. The total hux and PD variations seem 
to be uncorrelated, owing to the significant contributions of the unpolarized component to 
the total emission. However, our study shows the co-spatiality of emission at high and low 
energies. The observability of quasi-simultaneous emission provides hints about the location 
of emission region in the jet. The color variations and BWB trend indicate that it may likely 
be of shock and/or shock-initiated reconnection origin. The HMFM model, adopted to 
simulate the part of outburst, suggests that the magnetic reconnections more likely played a 
very important role in this event. The rising part of PA rotation is very well htted with “step 
phase” prohle. The same prohle also presents a reasonably good ht for PD variations, except 
for a slight deviation towards the end. This may be explained by including the contributions 
of other unpolarized emission components in the optical band. The excess in the ratio of the 
toroidal component to the poloidal component of the magnetic fields during quiescent and 
flaring episode indicate the helical component to contribute in the flux increment (Table [2]) . 
In conclusion, our study suggests that this outburst event is more likely to be the similar 





















to 2008 outburst state of S5 0716+714 , and is probably triggered by the shock-initiated 
magnetic reconnections taking place in the emission region in the jet. 
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Fig. 1.— Multi-wavelength lightcurve of S5 0716+714 showing the recent outburst activity 
during 2015, January. Figures 1(a) & 1(b) respectively represent the Fermi (>0.1 GeV) 
and X-rays (0.3 10.0 keV) light-curves, whereas the Figs. 1(c) & 1(d) present UV/optical 
magnitudes from Swift-UVOT and MIRO, respectively. Last two panels (1(e) & 1(f)) are 
PD and PA variations. The ‘Stew R’ notation in Fig. 1(d) stands for R band data from 
Steward Observatory.. 
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Fig. 2.— Left: R band and I,V, B optical bands correlations curves for fluxes in different 
optical bands. The X-axis shows the R band magnitude whereas the Y-axis represents 
magnitudes in other optical bands namely, I, B, & V. Right: The flux-color plot for optical 
observations. The X & Y axes are the R-band magnitude and (B-V) color, respectively. 
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Fig. 3.— Left: The PD variation with individual segments fitted by straight lines (LI to 
L5). The fitting is performed using least square fitting algorithm. Right: The PA variations 
fitted with an exponential rising and falling profile given by [f(x) j\ e ~a( a:—z0)j_ 

















18 


t 



3 t 



LOS 


Fig. 4.— Left: Sketch of the interaction between the emission region and the disturbance 
in the comoving frame of the emission region, at different epochs. The emission region is 
pervaded by a helical magnetic held and a turbulent component (only the helical component 
is sketched). The disturbance is stationary in the observer’s frame, but in the comoving frame 
of the emission region, the disturbance is then moving up with Lorentz factor T. The orange, 
red, green and blue regions refer to the locations of the disturbance before the hare (to), the 
rising phase (fi), peak (t 2 and £3) and declining phase (£4), respectively. Right: The red, 
green and blue shapes indicate the shape and location of the haring region, corresponding 
to the disturbance at t x ~ f 3 , respectively, observed simultaneously, taking into account the 
LTTEs. Since Z > 2 R, the peak state will stay for few hours or few days, depending upon 
the Z/R ratio ( Zhang et ah 2015). 
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Fig. 5.— Left: The model reproduced PD and PA, generated for the duration of MJD 57045.5 
to 57047.5. The PD variations overplotted with the modeled curve (thick black line). The 
bottom panel represents the data and model for PA swing. Right: The broadband SED of S5 
0716+714 for aforementioned duration. The multi-wavelength light-curves and respective 
models are not shown here as the flux remains almost constant (though highest) for the 
duration under consideration. 


Table 1: Swift XRT Spectral Fitting Parameters. 


ObsID 

Exposure 

(s) 

Time 

(MJD) 

Sp. Index (T ) 

log W (Fo.3-io.Okev) 
(erg cm ~ 2 s^ 1 ) 

x 2 

dof 

00035009145 

1043 

57023.21 

2.3 ±0.4 

-11.4 ±0.10 

4.6 

3 

00035009146 

821 

57029.01 

1.9 ±0.4 

-11.3 ±0.13 

1.6 

2 

00035009147 

976 

57041.09 

2.7 ±0.3 

-11.0 ±0.06 

5.0 

11 

00035009148 

1091 

57042.75 

2.4 ±0.2 

-10.9 ±0.06 

14.1 

13 

00035009149 

961 

57043.41 

1.8 ±0.6 

-10.7 ±0.20 

1.7 

3 

00035009152 

1348 

57044.02 

2.3 ±0.2 

-10.7 ±0.05 

23.1 

19 

00035009153 

6895 

57044.29 

2.3 ±0.1 

-10.8 ±0.03 

38.4 

50 

00035009154 

978 

57045.01 

2.5 ±0.3 

-10.6 ±0.06 

6.6 

9 

00035009156 

9574 

57045.14 

2.5 ±0.1 

-10.7 ±0.02 

92.2 

89 

00035009157 

1688 

57047.14 

2.6 ±0.1 

-10.6 ±0.03 

17.2 

26 

00035009158 

6557 

57047.22 

2.6 ±0.1 

-10.6 ±0.02 

79.8 

68 
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Table 2: The best estimates of the model parameters. 


Paramters 

Value (CGS) 

Bulk Lorentz factor 

20 

Length of the emission region Z (cm) 

6.06 xlO 16 

Radius of the emission region R (cm) 

2.25 xlO 16 

Length of the disturbance L (cm) 

6.06 xlO 15 

Radius of the disturbance A (cm) 

2.25 xlO 16 

Orientation of LOS (deg) 

90 

Electron acceleration time-scale (Z/c) 

5.50 xlO -3 

Electron escape time-scale (Z/c) 

6.00 xlO” 4 

Electron density (cm^ 3 ) 

21.7 

Helical magnetic held strength (G) 

0.5 

Helical pitch angle (deg) 

47 

Helical pitch angle during hare (deg) 

75.5 
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